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Abstract 

Among 12672 Feynman diagrams contributing to the electron anomalous magnetic moment at 
the tenth order, 6354 are the diagrams having no lepton loops, i.e., those of quenched type. Because 
the renormalization structure of these diagrams is very complicated, some automation scheme is 
inevitable to calculate them. We developed an algorithm to write down FORTRAN programs for 
numerical evaluation of these diagrams, where the necessary counterterms to subtract out ultra¬ 
violet subdivergence are generated according to Zimmermann’s forest formula. Thus far we have 
evaluated crudely integrals of 2232 tenth-order vertex diagrams which require vertex renormaliza¬ 
tion only. Remaining 4122 diagrams, which have ultraviolet-divergent self-energy subdiagrams and 
infrared-divergent subdiagrams, are being evaluated by giving small mass A to photons to control 
the infrared problem. 

PACS numbers: 13.40.Em, 14.60.Cd, 14.70.Bh, 11.15.Bt, 12.20.Ds 
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I. INTRODUCTION 


The anomalous magnetic moment of the electron, also called the electron g — 2, is one 


has been measured with steadily increasing precision 
electron and the positron available in the literature 


of the most fundamental quantities of particle physics_. Since its discovery in 1947 [1] it 


. The best values of g — 2 of the 


a e - (exp) = 1 159 652 188.4 (4.3) x 1(T 12 
a e + (exp) = 1 159 652 187.9 (4.3) x 1(T 12 


( 1 ) 


were obtained by the Penning trap experiment. Here a e = ^{g — 2) and the numerals in 
each parenthesis represent uncertainty in the last few digits of the respective values. The 
consistency of a e - with a e + in (D within the experimental accuracy exhibits that CPT is a 
very good symmetry of the universe. 

At present a new experiment is being carried out by a Harvard group using a new trap 


with cylindrical cavity 


|, which is capa 


with the help of analytical calculation 


dc of controlling the electron-cavity-wall resonance 
This experiment will reduce the measurement 
uncertainty of m substantially. It will enable us to test the validity of QED to a very 
high degree and to determine the fine structure constant a to an unprecedented precision 
of 7 x 10 10 or better, which is an order of magnitude better than the best non-QED value 
available at present pj. 

Of course, such a feat requires availability of theoretical calculation of matching preci¬ 
sion. Within the framework of the Standard Model the QCD and weak interaction parts 
of the corrections to a e are known to be so small that their uncertainties do not affect the 
determination of a even with the expected precision of the new Harvard experiment. The 
uncertainties due to the QED contributions induced by the virtual propagation of muon 
and tan-lepton beginning at the fourth-order (a 2 ) are also known to be negligible within 
the current precision. Thus the electron g — 2 within the experimental precision of our cur¬ 
rent interest is determined almost entirely by the electron-photon interaction and can be 
regarded as a function of a alone. 

The latest evaluation of a e in the Standard Model, including the hadronic vacuum polar¬ 
ization contribution, hadronic light-by-light-scattering contribution, the electroweak effect, 
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and small QED contribution from virtual muon and tau-lepton loops 


is 0 ] 


a e = 1 159 652 175.86 (0.10) (0.26) (8.48) x 10” 12 , (2) 

where the uncertainties stem from (i) the remaining numerical uncertainty of the ad-term []J, 
(ii) the crudely estimated uncertainty of the a 5 -term js|, and (iii) that of the best non-QED 
a available at present, which is measured by the atom interferometry combined with the 
cesium 1)\ line measurement by the frequency comb technique ji| , 


a-\h/M Cs ) = 137.036 000 3 (10) [7.4 ppb] 


( 3 ) 


An important byproduct of the study of the electron <7 — 2 is that a more precise a can 
be obtained by combining the measurement (JTJ) and the theory of a e , which yields 


a~ 1 (a e ) = 137.035 998 834 (12) (31) (502) [3.7 ppb], 


( 4 ) 


where the uncertainties 12 and 31 are due to the a 4 and a 5 terms, and 502 comes from the 
experiment ©. The new Harvard experiment of a e is expected to reach a precision an order 
of magnitude better than that of ©■ The a 5 term will then become the largest source of 
unresolved systematic errors for obtaining a from a e . Thus an explicit evaluation of the a 5 
term is urgently needed for further improvement of a(a e ). 

The pure QED contribution can be written as 


ae( QED) = 4 2 '0 + 4 4) 0 2 + 4 8, 0 s + 4 8, ( 


a\ 4 
n) + 


( 5 ) 


The coefficients are evaluated by the perturbation theory. By now the first four of them 
have been obtained ^j]: 

0.5, 


4 2) = 


4 4) = -0.328 478 965 • • • , 
Af ] = 1.181 241 456- •• , 

Hj 8) = -1.728 3 (35). 


( 6 ) 


4 10) has n °t yet been evaluated. An educated guess is that it may be found within the 
range (—3.8, 3.8) 

The first theoretical calculation of a e was carried out analytically by Schwinger in 1948 
jinj |. The number of Feynman diagrams involved was just one in this case. The excellent 
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agreement of his calculation with the measurement 1] was one of the pivotal triumphs of 
renormalization theory of QED, which was just being developed. Refinement of theory to 
the fourth-order involves seven Feynman dia gram s. It took more than 7 years before the 


analytic value of A[ 4> was obtained in 1957 


Analytic evaluation of A® is far more 


challenging requiring evaluation of 72 Feynman diagrams. It took effort of many physicists 
and many years of hard work and was completed only in 1996 jl3 | . 

The numerical evaluation scheme was developed by one of the authors (T. K.) and Cvi- 


tanovic for the evaluation of the sixth-order contribution 


to the eighth-order 


0 


0 Q 15] 


and extended later 


(o\ 

Up to now Aj ; was calculated only numerically, which involves 


the evaluation of 891 Feynman diagrams. Although the initial crude result was published in 


1974 [3], improvement of numerical precision re qui red many years of extensive computation 


|7]. From the viewpoint of obtaining A^ 10 ' 1 


and the final result was published only recently 
the numerical integration approach is the only practical choice at present. 

The contribution to the a 5 term of the electron g — 2 comes from 12672 vertex-type 
Feynman diagrams, which can be categorized into 6 sets according to their structures and 


classified further into 32 gauge-invariant subsets. See Appendix El and Ref. 


• 0 


for the 


details of classification. Most subsets that contain closed lepton loops have relatively simple 


. 0 


and 


structure and are calculable by a slight extension of the method developed in Refs. 
flf| . As far as the muon g—2 is concerned, they cover the subsets that give rise to the leading 
contributions. Thus far 17 of 32 subsets have been evaluated and reported in Refs. [3] and 


0 - 


10 . 


Detail of these works is presented in Ref. 

For the electron g — 2, however, none of 32 subsets is dominant so that all must be 
evaluated. A particularly difficult one is Set V, a huge set consisting of 6354 vertex diagrams, 
all of which have pure radiative corrections and no closed lepton loops (see Fig. EH in 
Appendix El)- Throughout this paper these diagrams will be referred to as “quenched-type 
(q-type)” since they are analogous to the so-called quenched diagrams of QCD. The difficulty 
of Set V stems from the fact that many of them have very large number of ultraviolet (UV) 
and infrared (IR) divergences. This makes the previous approach highly impractical since 
it runs into an extremely sever logistic problem. Unless this is solved, it will be close to 
impossible for mortals to deal with Set V and hence the entire tenth-order electron g — 2 
without making mistake. 

The purpose of this paper is to present our solution to this very difficult problem. We have 
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developed a scheme of automatic code generation which enables us to generate renormalized 
integrals for all diagrams of Set V with a breathtaking speed. Outputs of this code are ready 
to be integrated by numerical means. 

We begin by exploiting the equation derived from the Ward-Takahashi identity, which 
relates the sum of a set of vertex diagrams to a self-energy-like diagram Q|. This relation 
together with the time-reversal invariance of QED enables us to reduce the number of inde¬ 
pendent integrals of Set V diagrams drastically from 6354 to 389. The method starting from 
the Ward-Takahashi (W-T) summed diagrams is called Version A, while the conventional 

n 

approach by vertex diagrams is referred to as Version B for the sake of distinction [20|. 

The systematic scheme for constructing numerical integration code consists of the follow¬ 
ing steps. 


(I) Identify the diagrams contributing to the electron g —2 and their UV- and IR-divergent 
sub diagrams. 

(II) Carry out momentum space integration exactly using a home-made integration table 
and convert it into an integral over Feynman parameters. The result is expressed 
symbolically as a function of quantities U, By, and Aj, which are homogeneous poly¬ 
nomials of Feynman parameters. We call them building blocks. 


(Ill) Find the explicit forms of U and By which are determined from the topological struc¬ 
ture of the Feynman diagram Q obtained by removing all external lines and disre¬ 
garding the distinction between an electron line and a photon line. 


(IV) To prepare for numerical integration on computer UV and IR divergences must be 


removed from the integrand beforehand. In Ref. Q] a regularization scheme was 
developed in which divergences are eliminated by point-by-point subtraction by coun¬ 
terterms which are derived from the original integrand by simple power-counting rules. 
This scheme is denoted as //-operation for the UV divergence and /-operation for the 
IR divergence Q Q . 


(V) Counterterms thus constructed can be identified with only the UV divergent parts of 
the renormalization constants so that the result of step (IV) is not fully equivalent 
to the standard on-shell renormalization. The difference between full renormalization 
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and intermediate renormalization must therefore be evaluated by summing up all 
subtraction terms. 


The scheme (I)-(V) itself is completely general and applicable to any order of perturba¬ 
tion. Step (II) was quite difficult already in the sixth-order and still harder in the eighth- 
order. It was carried out entirely by computer with the help of algebraic manipulation 


programs such as SCHOONSCHIP |22| and FORM |23|. Steps (I), (III), (IV) and (V) were 
simple enough in the sixth-order case and still manageable in the eighth-order case to be 
executed by hand manipulation. 

It is evident, however, that such a pedestrian approach is no longer adequate for the 
calculation of the tenth-order diagrams so that a highly automated approach is required. 
This applies not only to the step (II) but to all other steps. It turns out that the automation 
scheme can be formulated quite efficiently for the q-type diagrams by making full use of their 
inherent properties. For Step (I) a systematic procedure for the generation of diagrams and 
the identification of UV divergent subdiagrams is possible. The graph-theoretical notions 
are easily identified for this type of diagrams, which enables the automated construction of 
topological quantities in Step (III). 

The UV subtractions in Step (IV) can be organized by following the Zimmermann’s forests 
of subdiagrams exactly The forests are constructed as combinations of subdiagrams 
identified in Step (I). 

This enables us to write a code which controls all steps (I), (II), (III), and (IV) auto¬ 
matically. Namely, we have obtained a code which turns an input of single-line information 
characterizing the structure of a Feynman diagram into a fully renormalized Feynman para¬ 
metric integral. 

Thus far we have obtained FORTRAN codes of renormalized integrals for 2232 vertex 
diagrams which contain vertex renormalization subdiagrams only. Crude evaluation by the 
Monte-Carlo integration routine VEGAS shows that our scheme works well as expected. 

The next step is to evaluate the remaining 4122 diagrams which have self-energy subdi¬ 
agrams. These diagrams have also IR divergences. The simplest way to deal with the IR 
problem is to give a small mass A to photons, which can be implemented with a minor ex¬ 
tension of the automating code. Of course, the numerical result will have an uncertainty of 
order A. It may also suffer from non-negligible digit-deficiency problem commonly encoun¬ 


tered in numerical integration 


Nevertheless it will be good enough for getting crude 


6 




values of A^ 101 so that we follow this approach as the first step. 

This scheme has thus far been tested successfully with the sixth-order q-type diagrams 
and reproduced the analytic result after proper treatment of residual renormalization terms. 
The next step is to check the eighth-order q-type diagrams, which is numerically known. It 
seems successful so far except for a few diagrams which have severe IR divergences. Those 
diagrams may require minor modifications to the present automation code. After these 
exercises we will tackle the tenth-order problem. 

To obtain a result independent of A it is necessary to extend the automating code to 
include IR subtraction terms constructed in a manner similar to that of UV counterterms. 
To complete this calculation we must evaluate the contribution of residual renormalization 
terms, which consists of integrals of up to eighth-order for 6804 UV-divergent subdiagrams. 
The result of these works will be reported in the subsequent papers. 

This paper is organized as follows. In Section ITT1 we briefly review the parametric integral 
formalism to obtain the anomalous magnetic moment of leptons in Version A approach. 
In Section nm we describe A-operation for the subtraction of UV divergence which derives 
from a single UV-divergent subdiagram. We then proceed to the case involving more than 
one subdiagram in Section II VI in relation to the forest structures. It is one of the essential 
ingredients for our automated scheme. In Section 0 and EU we focus on a specific type of 
diagrams, namely, q-type diagrams. We discuss their intrinsic properties in Section M We 
see in Section EU that use of these properties allows us to develop algorithms to identify 
UV divergences in terms of all relevant forests. It enables us to construct the renormalized 
amplitude in an automated manner. In Section lUlIl we show the whole flow of our automated 
scheme in detail. Section IVUII is devoted to the conclusion and discussions. In Appendix 1X1 
we show the classification of diagrams that contribute to A^ . Appendix El contains useful 
formulae for computing the basic building blocks of the Feynman integrals which can be 
readily adapted to the programming languages such as C++ and FORTRAN. 


II. GENERAL FORMALISM 


In this section we present the general formalism for evaluating QED contribution of 
anomalous magnetic moment of lepton in perturbation theory. It is a brief summary of the 
literature [Ijj |2l|, which is included here to provide concrete prescription of formulation 
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for our automation scheme presented in the latter part of this paper, and also to make this 
article self-contained. 

A. Anomalous magnetic moment of lepton 

The magnetic property of a lepton can be studied through examining its scattering by a 
static magnetic field. The amplitude of this process including interactions with the virtual 
photon fields can be represented as follows, by taking account of the gauge symmetry, 
invariances under Lorentz, C, P, and T transformations: 

eu(p") + ^-a^q u F 2 (q 2 ) u(p') A*(q) , (7) 

where p' = p — q/2, p" = p + q/2, q = p" — p' and (T ,JV = ^( 7 ^ 7 ^ — 7^7^)- A* is the vector 
potential of the external static magnetic field. F\ and F 2 are called the charge and magnetic 
form factors, respectively, and the charge form factor is normalized so that i*i( 0 ) = 1 . 

The anomalous magnetic moment a e is the static limit of the magnetic form factor F 2 (q 2 ), 
and it is expressed as 

a e = 7 * 2 ( 0 ) = Z 2 M (8) 

with 

M = lim TrfP^p, q) T") , (9) 

P —>0 

where Z 2 is the wave function renormalization constant, Y v is the proper vertex part, and 
P„(p,q ) is the magnetic projection operator, 

q) = jp-j + rnj rrii'p 1 - ^m 2 + p v ^ + rnj . (10) 

Here, the momentum of incoming lepton p — \q and that of outgoing lepton p + \q are on 
the mass shell so that p and q satisfy p 2 = m 2 — \q 2 and p ■ q — 0. 

We evaluate the anomalous magnetic moment a e in the framework of perturbation theory. 
Because of renormalizability of QED it can be written as a power series in - whose coefficients 
are finite and calculable quantities. 

B. Construction of Feynman parametric integral 


In perturbative analysis of QED the amplitude is usually expressed as an integral of loop 
momenta flowing through the Feynman diagram. In this paper we convert it into an integral 



of Feynman parameters z % assigned to internal lines |13j, |26| • 

We consider a 2nth-order lepton vertex diagram Q which describes the scattering of an 
incoming lepton with momentum p — q/2 into an outgoing lepton with momentum p + q/2 
by an external magnetic field. Q consists of 2n+l interaction vertices connected by 2 n lepton 
propagators and n photon propagators, which are given in the form (in Feynman gauge): 

( 11 ) 


+ rrii —ig^ u 


2 9 J 

Pi ~ mf 


2 2 ? 
Pi ~ m i 


respectively. The momentum p,- may be decomposed as p* = &* + <p, in which fc* is a linear 
combination of loop momenta, while is a linear combination of external momenta, m; is 
the mass associated with the line i, which is temporarily distinguished from each other. 

We introduce an operator D % /J by 


Dr = - 


fCXD Q 

dm 2 


dq, 


( 12 ) 


z/x 


and replace each numerator ]fi — fa + fa of lepton propagators CP by Ipi. Since DP does 
not depend on k t explicitly, the numerators can be pulled out in front of the momentum 
integration as far as the integrand is adequately regularized. 

The product of denominators are combined into one using the formula, 


N 


1 


i=1 *- 1 


' N i 

II/ dZi 

i= 1 J o 


N 


di-E 


i =1 


N 


N 


(13) 


Y^ z,Xi 


, i= 1 


The sum JT ZiXi is a quadratic form of loop momenta so that it can be integrated analyti¬ 
cally. As a consequence the amplitude is converted into an integral over Feynman parameters 
Zi which is expressed in a concise form as 

1 


r 5 = (!) (n - 1)! F "j ( dz) g 


U 2 V r ‘ 


where N = 3n and 


N 


N 


( dz)g = dZi 5 I 1 - ^2 Z i ) > 

i=1 \ i=l J 

N 

v = Zi ( m i 2 - Qi ■ Qi) 1 


1=1 


N 


«!' = - 5 E 


3 = 1 


B 'ij = B i) - y 


U 


(14) 

(15) 

(16) 

(17) 

(18) 
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In Eq. ED we have omitted the factor ( a/n) n for simplicity. U and Bij are homogeneous 
polynomials of degree n and n — 1 in Feynman parameters {zi}, respectively. Their precise 
definitions are given in later sections. The operator F 1 ' is of the form 

n 

F" = 7 Q1 C0i + m ib“ 2 • • • 7 " • ■ ■ 7“ 2n - 1 ijp2n + m 2n ) 7 Q2 ™ J] , (19) 

k =1 


where TT, o a . a . is a diagram-specific product. If Q has closed lepton loops F 1 ' also contains 
appropriate trace operations. 

Note that can now be brought into the z- integral. The operator D i ,l in F y acts on 
l/V n as 


1 0' p 
jj 


V r 


V r 


i Q'.^Q'y 

n v 1 = 

i j yn y 


g ,,v B{ 




DfD "D k p -i- = 

1 3 K \Tn 


2 (n - 1 ) 


2 (n - 1) UV "- 1 ’ 

QTQ'j'Q'k 

yn 

{gTB^M + 9 up B' k Q'S + g p »B' ki Q'n 


UV 


n —1 


( 20 ) 

( 21 ) 


( 22 ) 


The result of this operation may be summarized as a set of rules for a string of operators 



a) when Ip t and Ipj are “contracted”, they are turned into a pair of 7 ^ and 7 ^ times a 
factor (—| B^). 

b) uncontracted D, is replaced by Q[. 


As a consequence the action of F" produces a series of terms of the form 


F" 


U 2 V r - 


+ 


F? 


U 2 V n U 3 V n ~ 1 


+ 


(23) 


where F k are polynomials of B[, and Q\. The subscript k denotes the number of contractions. 


F k also includes an overall factor 


(n — l)(n — 2 ) • ■ ■ (n — k ) 

It is convenient to replace vectors Q p by scalar functions. Suppose the momentum 


p^ —— enters the graph Q at the point A, follows the path V = V(AC), and leaves at C, 
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FIG. 1: A vertex diagram with given external momenta. A choice of paths V' and V" are shown. 


o p 

and p^ 1 + — enters at C, follows the path V" = V(CB), and leaves at A. (Fig. (JTJ).) This 
can be expressed concisely by 


<?/ = Vjv 




+ VjP" 



(24) 


where ippi = (1, —1,0) according to whether the line j lies (along, against, outside of) the 
path V'. Similarly for pj-pn. Substituting Eq. EH in Eq. (fT7l) we obtain 

QT = Ar (p“ - A ) + Af (p- + q -Aj . (25) 

where 


N 


AT' = 


u 


J2Vj-P'ZjB'ji■ 


(26) 


3 = 1 


qb 


qb 


Similarly for AT ■ AT , AT will be called scalar currents associated with + -E-, 

respectively. 

If we choose a path V = V(AB) for the corresponding scalar current becomes AT = 
AT' + AT" ■ Note that the choice of V(AB) is flexible as far as the end points A, B are 
fixed. Note also that V(AB) no longer depends on C. 


C. Building blocks, B tJ and U 

In our formalism, the parametric functions B a p and U provide the basic building blocks 
which are defined on the chain diagram corresponding to the diagram Q. Here a, f3 refer 
to the chains; a chain is a set of internal lines that carry the same loop momentum. The 
chain diagram is derived from Q by amputating all the external lines and disregarding the 
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1+a 




FIG. 2: A diagram (left) and the chain diagram derived from it (right). 

distinction between the types of lines. Every chain is assumed to be properly directed. B a p 
and U are homogeneous polynomials of degree n — 1 and n, respectively. They are the 
quantities that reflect the topological structure underlying the diagram Q. 

B a p and U can be obtained recursively by the following relations, 

B a /3 = ^2 £/ 3 ,c Ug/c , (27) 

C 

i\, s U = ^2 £\,s z aB\ a , for any A e s , (28) 

a 

starting from U = a for a single loop. Here the summation over c runs over all self- 
nonintersecting closed loops on Q. The loop matrix £ Q:C is a projector of chain a to the loop 
c, which takes (1, —1, 0) according to whether a is (along, against, outside of) c. Ug/ C is the 
U function for the reduced diagram Q/c that is obtained from Q by shrinking the loop c to 
a point. The loop s in Eq. (j2Hj) is an arbitrary closed loop. 

Alternate and equivalent formulae for B a p and U are obtained in the following manner. 
Suppose a set of independent self-nonintersecting loops (called a fundamental set of circuits) 
is given and define U s t by the summation over all chains by 

U s t ^ £,a,s £ a,t : (29) 

a 

where s,t are labels of circuits in the set. Then, U and B a p are given by 

U = det U st , (30) 

st 

B a p = uY J Us^t(U- 1 )st. (31) 

st 

For a given diagram Q, first we have to identify the fundamental set of circuits, and construct 
the loop matrix £ a<s . Then we can obtain U and B a g according to the formulae above. 


12 



Bij of the lines i , j is identical with B a p whose indices are such that i G a and j G f3. 
Bij satishes a so-called junction law on each vertex if the diagram Q were regarded as an 
electric circuit in which the Feynman parameter z t corresponds to the resistance of the line 
v. 

= ° (32) 

i 

for any vertex v and any internal line j, where e m is called incident matrix dehned by 


1 if the line i enters the vertex v, 

< —1 if the line i leaves the vertex v, 
0 otherwise. 


(33) 


Bij also satishes a loop-law given by the following relation for arbitrary closed loop s and 
arbitrary line j: 

= 0. (34) 

i 

These relations reduce the number of independent elements among Bij. It also provides 
consistency checks which are useful in the actual calculations. 


D. a e from a set of vertex diagrams summed by Ward-Takahashi identity 


A set of vertex diagrams which are derived from a self-energy diagram by inserting an 
external vertex in every lepton propagators share many properties. Actually we can even go 
further to relate those integrals to a single integral of the self-energy-like diagram through the 
Ward-Takahashi identity. This relation is useful when we consider higher order calculations 
because it reduces the number of independent integrals substantially. 

It is well known that the proper vertex T M = + A M and self-energy part E are related 

by the Ward-Takahashi identity 




= -E [p + -q + E [p - -q 


(35) 


This relation holds perturbatively as well for Eg representing the lepton self-energy diagram 
Q and the sum of vertex diagrams Ag that are obtained by inserting an external vertex into 
Q in every possible way. Differentiating both sides of Eq. © with respect to and taking 
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the static limit q —> 0 of the external magnetic field, we have 


A v ip, q) ^ -q 1 * 


dA»(p,q) 

dq v 


Jq=0 


gS(p) 

dp u 


(36) 


We may evaluate a e starting from either side of this expression; a straightforward way is to 
calculate each vertex diagram individually and to gather them up according to the left-hand 


0 ), 


side (Version B approach in Ref. l20||), or else we can combine the set of vertices into one 
according to the right-hand side (Version A approach). We adopt Version A in the present 
study. 

In the Feynman parametric form, the 2nth-order magnetic moment associated with a 


self-energy-like diagram Q can be written as 


M<2 " ) = T r n - 1)1 {izh —Tpv^ + < n + z w 


e + c 


where E, C, N, and Z are a set of operators defined as 


N = -TV[Pfo,(2GF)] , 


1. 


E = jTr[P;'E„] , 
C = iTr[PfC„„] , 
Z = iTr[Pf Z„„] . 


(37) 

(38) 

(39) 

(40) 

(41) 


The magnetic projectors P” and P^ u are derived from Eq. (1101) by averaging over the direc¬ 
tion of and take the following forms: 


p r = V _( 1 + t£W 


p ^ _ 
- r 2 ~ 


1 + 


m 


3 m J m 


^ M v P n 

g M — 7^7 H-7-7^ 


m 


m 


(42) 

(43) 


The operator F is the numerator part of the self-energy-like diagram Q constructed in the 
similar form as Eq. PI): 


F = 7 ai (^i + ™ib“ 2 • • • 7“ 2 "' 1 (^2n-i + m 2n . 1 ) 1 a ^Y[g 0 


l fc u Ok 


(44) 


k =1 


which may contain appropriate trace operations if Q has closed lepton loops. The operator 
E" is defined by 


W = 


d¥ 

dPu 


E A ‘ F ‘“’ 

all leptons 


(45) 
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in which FT is obtained from F by substituting in the ?'th line: 


(9i + TOj) -* 7 ". 


( 46 ) 


The operator is dehned by 

^ = (47) 

3 

The sum runs only over the lepton lines into which the external photon line can be inserted. 
Z T 17 is obtained from F by substituting in the jth line: 

{JPi + m j) -»• \ [7 M 7 v (Jpj + m j) - (Jpj + • (48) 

The operator C ^ is dehned by 


C #w = J]cyF i f', (49) 

i<j 

where i and j refer to all lepton lines. C tJ is given by 

C« = 4 E , (60) 

k<l 

where k, l are taken from the lepton lines that belong to the path on which the momentum 
q u of the external magnetic held hows. FT*' is obtained from F by substituting in the ith 
and jth lepton lines: 

(JPi + rrii), {Jpj + mj) -> Y . (51) 


G is given by 

G = J2 z iAi, (52) 

l 

where the summation runs over the lepton lines on which the external momentum p M hows 
(depending on the choice of path V(AB) for the scalar currents). 

We can now construct the integrand in the following two steps. 


(I) Express the integrand as a function of symbols B tJ , A iy U, V, and 

(II) Express those building blocks explicitly in terms of the Feynman parameters Z{. 

Step (jH) can be achieved analytically by algebraic manipulation programs such as FORM 
23l | . All the integrals are generated from a small number of templates with the permutation 
of indices according to the specihc structure of each diagram. Step is performed along 
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the prescriptions outlined above, once and U are obtained by the formulae in Section lll Cl 


The magnetic moment contribution P7j) now can be expressed as a parametric integral: 



M {2n) 


(53) 


III. SUBTRACTIVE UV RENORMALIZATION PROCEDURE 

The amplitude thus far constructed in the previous section is divergent in general, and 
the divergences must be removed before carrying out the integration numerically. The UV 
divergence arises when one or more loop momenta go to infinity. This is seen in Feynman 
parameter space as the parameters z t that belong to loops in a subdiagram go to zero 
simultaneously. It allows power counting rules for identifying the emergence of singularities 
in a similar manner to the ordinary momentum integration. 

We adopt here the subtractive on-shell renormalization. In this scheme the renormaliza¬ 
tion term involving an mth-order vertex renormalization constant L m is given of the form 
—LmMn-m, where M n ^ is a g —2 term of order n—m. The renormalization constants that ap¬ 
pear in QED are the mass renormalization constant 5m, the wave-function renormalization 
constant B, and the vertex renormalization constant L. They are determined on the mass 
shell, and thus the coupling constant e and lepton’s mass m are guaranteed to be physical 


ones. 


To perform renormalization numerically our strategy is to prepare the subtraction term 
as an integral over the same domain of integration as the original unrenormalized amplitude, 
and to perform point-wise subtraction in which singularities of the original integrand are 
canceled point-by-point on the parameter space before the integration. To achieve this the 
renormalization constant L m and the lower-order g —2 term M n _ m are both expressed in the 
parametric integral and combined by the Feynman integral formula. It is found, however, 
that the integral is intractable if L m is treated as a whole. Instead, we adopt the following 
two-step intermediate renormalization, in which L m is split by 



( 54 ) 


and only the UV-divergent part is subtracted. 
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The subtraction term M n ^ m is found to have a term-by-term correspondence with 
the UV-divergent term of the original integral M n , and thus cancels the UV singularities. It 
is identified from the original integrand by simple power counting rules. This procedure is 
formulated as K -operation. The treatment of the UV divergence of self-energy subdiagram 
is slightly more complicated. See Ref. HQ and Eq. dUD for details. 

The UV-finite part of the renormalization constant is treated separately together with 
those from other diagrams 1 . This step is called the residual renormalization. 

In this section we shall describe how to construct the intermediate renormalization term 
via A"-operation. It is shown that the subtraction term factorizes exactly into the UV- 
divergent part of the mth-order renormalization constant and by construction. This 

feature is crucial for the subsequent operation when the UV divergence arises from more than 
one divergent subdiagrams. Such cases are treated more thoroughly in the next section in 
relation to the forest structures. The factorization property is also significant for the residual 
renormalization step in the sense that the highest order of the residual part decreases by 
two, e.g., for the tenth-order diagrams it is sufficient at most with the eighth-order terms. 
Therefore the evaluation of the residual part reduces to lower-order integrals. 


A. UV divergent subdiagram 


The UV divergence associated with the subdiagram S is caused by the simultaneous 
limits ki —> oo of all loop momenta /q, i G S. In the parametric representation (1531) this is 
translated into the vanishing of the denominator U at a boundary of Feynman parameter 
space where 2 


\0{e) 

Zi = < 

0 ( 1 ) 


i e S , 
otherwise, 


(55) 


with e —> 0. 

To find how a UV divergence arises from a subdiagram S consisting of internal lines 


1 It may contain IR divergences in general and they are also subtracted in a similar manner as UV diver¬ 
gences, though this subject is not covered here. In this article we instead introduce cut-off to treat IR 
problems. 

2 The overall divergence of a self-energy-like diagram drops automatically after projecting out the magnetic 
moment contribution. 
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and ns loops, consider the integration domain In the limit e —> 0, the homogeneous 

polynomials in the integrand behave as follows. (See Section fill Cl for proofs.) 


U = 0{e ns ), V = 0{ 1), 


(56) 


B,j = 


(57) 


and 

0(e ns ~ 1 ) if i,jeS, 

0(e ns ) otherwise . 

Let ms be the maximum number of contractions of operator D, within S. Simple power 
counting shows that the m-contracted term of M^ 2n) in Eq. (|531) is divergent if and only if 


Ns — 2ns < min (m, ms) , (58) 

where min (m,ms) means the lesser of m and ms- If S is a vertex part, we have Ns = -iris 
and ms = ns- If S is a self-energy part, we have N s = 3 n s — 1 and m$ = n$ — 1. In both 
cases Eq. dSHD is satisfied only for m > ms- Let us denote the UV limit m of U and Bij 
as [^]&V anc ^ [Bij]uv 


B. if-operation 

We are now ready to set up the rules of A"-operation for constructing the intermediate 
renormalization term. Firstly, we summarize our notation. QjS denotes a residual diagram 
which is obtained from Q by shrinking a subdiagram S to a point. Q — S denotes a diagram 
obtained from Q by eliminating all lines that belong to S. 

The K -operation Kg is defined as follows. 

(1) In Eq. (j53j) . collect all terms which are maximally contracted within the subdiagram S. 

(2) Replace U, B %] , C l3 , and Ai appearing in the integrand with their UV-limits, [R]uv ; 
[Bij] uv, [Cp]uv, and [A]uv, respectively. 

(3) Replace V with Vs+Vg/s, where Vs and Vg/s are V functions of S and G/S, respectively. 

(4) Attach an overall minus sign. 

A naive UV-limit gives V — > Vg/s instead of step (JHJ). Since Vs is a higher order term in e, 
its addition in step © does not affect the UV-limit. But it is crucial because it enables us to 
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FIG. 3: A closed loop c running in Q. 





satisfy the exact factorization of the renormalization constant and the rest of the amplitude 
required by the standard renormalization ju|. Furthermore, it enables us to avoid the 
spurious IR divergence which Vg/s alone might develop in other parts of the integration 
domain. 


C. UV-limit of building blocks U, and C \ 




Let us now describe step by step how the building blocks of the integrand behave in the 
UV-limit (1551) . It is found that each of them factorizes into two parts, one of which depends 
solely on the subdiagram S , and the other on the residual diagram Q/S alone. Since the 
description given in the literature is somewhat sketchy, we shall fill in the gaps here in 
preparation for automation of the procedure. 

The U function is a homogeneous polynomial of Feynman parameters of degree n defined 
by Eq. (USD, which has a simple behavior in the limit (1551) |27fl 


[U] s vv = u s u g/s (=o(€" s )). 


(59) 


In order to obtain the UV limit of B^, let us note that for i 6 a, j e (3, B a g of Eq. (1271) 
can be written as 

Bij — ^ €i,c£j,cUg/c • (60) 


Since (Q/c)nS = S/(cf)S) and ( Q/c)/S = C//(cU<S), the UV-limit of Ug/ C becomes 

[Ug/d uv = U s /(cns) Ug/(cus) ■ 


( 61 ) 


The explicit form depends on how the loop c runs in Q : 
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c is contained in S. (i.e. c C S. ) See Fig. |3] (jgj). 

In this case S/(cC\ S) = S/c and Q/(cU S) = Q/S. Therefore 

Ri/jSv = U s/c U g/s (=0(£"“- 1 )). (62) 

The power of e decreases by 1 since S/c has one less loops than S. 

c runs outside of S. (i.e. c C (Q — S). ) See Fig. |H| (fbj). 

In this case5/(cfl5) =S and^/(cU5) = (Q/S)/c. Therefore 

[CfyJuv = Us U (e/s)/ c (= 0(e" 5 )) . (63) 

c is contained in both S and Q — S. (i.e. c fl S ^ 0 and c fl (£? — S) ^ 0. ) See 
Fig.Udsj). 

In this case cfl5 is an open self-nonintersecting path within S. It does not change 
the number of loops in S when the path is shrunken to a point. Therefore the 
scaling behavior is 

Kl/Jov = 0(€"“), (64) 

though the exact factorization does not occur. 

From these observations and Eq. (EOD we find the following behavior of B tJ in the UV 
limit. 

I) Bij for i,j e S. 

The closed loops appearing in the sum in Eq. © fall into either of the cases (jsj) or 
n, the former gives the leading contribution whereas the latter does not in the limit 
Thus we have 

[-^h'luv 'y y ^i-c' Cj,c/ Us/c' Uq/s 

c'cs (65) 

= Bfj Ug/s , 

where the superscript S denotes that Bf- is the 5-function defined on the subdiagram 

5 . 

II) Bij for i,j e Q/S. 

The closed loops appearing in the sum in dSDD fall into either of the cases © or fcj), 


Case (a) 


Case (b) 


Case (c) 
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both of which give the same order of contributions: 

[-®h']uv = 'y 1 £<i,c' £j,c’ Ug/ C i + 'y ^ £,i,c" £,j,c"Ug/ c " . (66) 

d in case JbJ d' in case q 

In the first term on the right-hand side the sum over closed loops d C (Q — S) is 
equivalent to the sum over loops in Q/S — {s}, namely the loops in residual diagram 
Q/S that does not pass through the point s, where s denotes a point into which the 
subdiagram S has shrunk. Therefore, the first term becomes 

Us ^ £i,cdj,c'U(g/ S )/c' ■ (67) 

c’c(g/s- W) 

In the second term the closed loop c" passing through the points A, B G S fl (Q — S) 
is decomposed into two open paths V(AB) = c" fl S and V(AB) = c" fl (Q — S). The 
sum over c" becomes the sum over a choice of points A, B and open paths V(AB), 
V{AB). It is shown ( 2 ^] that U$/v satisfies 

U s = U s/v- (68) 

P(AB) 

On the other hand the path V'(AB) becomes a closed loop in Q/S that passes through 
the point s to which S has shrunk. Thus the second term becomes 

Us ^2 &,c" i].c" U(g/sy c " ■ (69) 

c"Cg/S, c"Bs 

From Eqs. (R771) and (RTHl the UV-limit of B l3 is 

[Bd w = Bf/ S U s i,je(G/S). (70) 


III) B m j for m G S and j G Q/S. 

This case is relevant only when S is a self-energy subdiagram, since for the vertex 
subdiagram case the leading contribution comes from the terms in which all lepton 
lines in S are contracted with each other. 


We denote the lines which are attached to the subdiagram S by i and i'. (See Fig. Q]) 
The closed loop c that contains the lines mGiS and j G Q/S passes through i and i!. 
The sum over loops c is decomposed into the sum over V = cDQ and V = cD (Q — S). 


It is shown 


& 


that 


J2UvU s/ v = U s A s m , ( 71 ) 

v 
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FIG. 4: A self-energy subdiagram S and a closed loop c that passes through m G S and j G Q/S. 



where A ^ is a scalar current on the line m of the diagram S. The path V turns into 


a closed path d after shrinking S to a point which passes through the line i G Q/S. 
Therefore B m j becomes 



(72) 



The UV limit of the scalar current Aj follows from Eq. (121 where the path V (which 
replaces P') is taken arbitrarily between two points attached to external lines. We can 
always choose the path to avoid the line j so that B'j in Eq. (El becomes B tl . 

When S is a vertex subdiagram, it is sufficient to consider only Aj with j G Q/S, since in 
the leading contributions of the integrand all the lepton lines in S are contracted and there 
is no operator Di left to be turned into scalar current. The sum in Eq. El consists of two 
parts, one from V = V fl S and the other from V" — V fl (Q — S). In the limit El the 
scaling behavior m shows that the former part gives sub-leading contribution. Therefore, 
using Eq. d, we obtain 



[Aj] 


(73) 


When S is a self-energy subdiagram, the scalar currents of both j G Q/S and j G S are 
relevant. 


Case (a) Aj for j G Q/S. 


The same argument of the UV limit as in the vertex subdiagram applies to this 
case, which leads to 
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Case (b) A m for m E S. 

We choose the path V so that it avoids S. Then all B im in 
cm whose UV limits are given by Eq. (72). Therefore, 


[A. 


1 


m Juv 


U 


E 


= A 


G / S k&V 

g/s A s 


„ TD S/s a s 
Uik ^m 


fall into the type 


(75) 


where i is the line adjacent to S. 


We recall that CV,- is derived from the part 

(76) 

9=0 

of Eq. (1361) with the external vertex inserted into the line j and differentiated with respect 
to the external momentum q v flowing through the line i. When S is a vertex subdiagram, 
Cij for i or j in S have no overall UV divergence, since S has, effectively speaking, four 
legs: one photon line attached to the external vertex, the other internal photon line that 
is connected to G/S, and two internal lepton lines. So it is sufficient to consider the cases 
i,j G G/S, in which the UV-limit of C l3 becomes 

lcy " v= tW cr ' (77) 

When S is a self-energy subdiagram, the definition 0 of Ci 3 and the UV-limits of -By- 
lead to the following forms. 



[Cu-W - Ug/s c jk 


G/S 


[Cfg] uv — jj s Cfg + 


U S 


4 5 E^a-^/E^) 

k hCiS hCiS J 


j,k e G/S , 

1 


E z iK 


'G/S 


(78) 

f19 £ S, 

(79) 


[Cf. 


^ i<S /~i G/S 1 


jJuv 


U 


G/S 




g£S 


f9 Uc 


G/S 


>G/S 
J jk 


f e S,j e G/S, (80) 


k£Q/S 


where i is the line adjacent to S. 


D. Factorization property of UV subtraction term 

Now we proceed to examine the UV subtraction term along the steps of 77-operation to 
see that it factorizes into two parts. For simplicity we consider a vertex part Tg defined in 
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Eqs. m and (I2H1) . though the arguments apply to the general cases. 

Suppose the UV divergent subdiagram S is a vertex subdiagram. In step (1) of K- 

operation we pick up the terms which are maximally contracted within S. Such a term 

F k 

among the terms with k contractions, _ , have the form: 


JJ2+kyn-k ' 


U 2 V 


n—k 


B 


u 

i,jes 


U 




(81) 


V 


The first factor in the braces is a product of B^s with i, j G S , while the second factor is 
a product that consists of ( k — ns) Bi^ks and several scalar currents whose indices i',f are 
in Q/S. 

In step (2) we consider the UV limit (I55|h It is achieved by replacing the building blocks 
U, Bij and Aj by their UV limits, [I/]u V , [B tJ ]f :y , and [A,-]u V , respectively. Then Eq. (1HT1) 
turns into 


1 

7/2 

u s 


B s 
U s 


r>G/S' 


TJ 2 
u g/s 


U 


■■■A 


Q/S 


Q/S 


V, 


n—k 

Q/S 


(82) 


=s[<5] =g[Q/S\ 

The first part depends only on Zi with j 6 5, which we denote by g[S]. The second part 
depends only on with % G Q/S. It is denoted similarly by g[Q/S\. In the naive UV limit 
V leads to Vq/s- 

In step (3) Vg/s is replaced by Us + Vg/s- The integral now becomes 

1 


{dz) g g[S\g[G/S\ 


, n—k 


{Vs + Vg/s) 

We shall see that it factorizes into S and Q/S parts. Firstly, the identity 

" 1 dt 


(83) 


i — l-« 

./n S V S 


(-¥)• 


(84) 


is inserted into the integral, where zs and zg/s are defined by zs = Zi and zg/s = Zi, 

i£S i€iQ / S 

respectively. Secondly, we rescale the Feynman parameters as follows: 


Zi —> sZi i G S 

Zi —* tzi i G Q/S 


(85) 


Since V-functions are homogeneous polynomial of degree 1, they scale in such a manner as 
Vs —> s Vs and Vg/s —■ y tVg/s■ Other parts of the integrand and the integration measure 
also scale accordingly. 
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Then using the Feynman integral formula: 


T (k + l) / dsdt5(l — s—t) 


s k 1 t l 1 


_ rW r(Q 

(.sA + tB)*+' A k B L ’ 


( 86 ) 


the integral is shown to be factorized into two parts: 


dz s 5(l-z s ) g\S] / dzg/sS(l-zg/s) g[G/S] / dsdt5(l-s-t) 


s a-l t P-l 


(sV S + tVg/s) 


a+/3 


= f {dz)s x [ ldz k/s 


3 \SJS\ 

V 13 

V Q/S 


(87) 


where a and (3 are constants determined by the rescaling (USD- 

Based on those observations the whole integral of the vertex part Tg is shown to be 
factorized in UV limit as 

K 5 r/ = L^rg/s , (88) 


where L^ IV is the UV divergent part of the vertex renormalization constant L$ and r g/ S is 
the vertex part of the residual diagram G/S. 

When S is a self-energy subdiagram, the factorization is not apparent because not all 
Ip m with 77i G S are contracted. From Eqs. (1721) and m we can symbolically write the 
uncontracted Jp m as 


[WSv = ^# /S , (89) 

where i" is a fictitious line related to i and . After a little algebra one finds 171 171 


K s r e - = dm" v r e y (i ., 


+ fi'.( v r 


g/s,i' 


(90) 


where is the UV divergent part of the mass renormalization constant 5ms and 

is that of the wave function renormalization constant Bs■ G/S{i*) denotes the diagram 
obtained by shrinking S to a point, where i* indicates two-point vertex between lines i and 
i!. G/S,i' denotes the diagram derived from G by shrinking S to a point and eliminating 
the line i'. It can be reduced to the form TJU after integration by part with respect to 
The factorization of if-operation is crucial when there are more than one subdiagrams 
that cause UV divergences, since this property guarantees that the successive operation of 
is consistent with the forest structure. 
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IV. FOREST FORMULA 


A Feynman diagram that appears at higher-order terms of perturbation theory may have 
complicated UV-divergence structures. In many textbooks they are treated in a recursive 
formulation so that the inner sub divergences of a renormalization part should be subtracted 
prior to the subtraction of its own overall divergence. It is natural in the framework of 
renormalization theory, for it is derived from the requirement that the divergences should 
be resolved by local counterterms. It is also tractable in general for hand manipulations 
since the subtractions are performed step by step from lower order parts and the number of 
steps are, as it turns out, relatively small. It is, however, not so convenient in our numerical 
approach in which the singularities due to the UV divergences are canceled point-by-point 
in the Feynman parameter space. To achieve this we have to prepare the subtraction terms 
as integrals defined in the same parameter space as that of the original unrenormalized 
amplitude. 

An explicit solution of the recursive formulation is given by Zimmermann’s forest formula 


24]. Each source of the UV-divergence is related to a forest , a set of UV-divergent subdia¬ 


grams, and the subtraction term associated with the forest is constructed by the subtraction 
operations for the subdiagrams applied successively to the unrenormalized amplitude. The 
whole subtraction terms are generated along the complete set of forests. 

In our numerical approach the subtraction operation is given by R-operation for a single 
sub divergence. As seen in the previous section it retains the factorization property, which 
guarantees the successive K -operations. Therefore, once a UV-divergent structure is known 
in the form of a forest, we can obtain the integrand of the subtraction term 15|. Although 
the subtraction scheme presented here is identical with that developed in Ref. [l51, it is more 
readily adaptable for code generation. 

The forest formula is much more useful for the automated scheme. The forests are given 
by the combinations of non-overlapping subdiagrams. So the complete identification of 
UV-divergent structures is obtained by purely combinatorial procedure from the set of all 
UV-divergent subdiagrams. Thus it is readily implemented in terms of forests, which enables 
us to obtain fully UV-renormalized amplitude of a diagram Q. 
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A. Definition of forests 

We begin with the inclusion relations between subdiagrams. When two subdiagrams S t 
and Sj share no vertices nor lines they are called disjoint When all lines of S t belong to 
Sj, the subdiagram S t is included in the subdiagram Sj. In this case S t and Sj are called 
nested. Otherwise they are called overlapping , in which Si and Sj share some vertices and 
lines while one is not included in the other. 

A forest is defined as a set of subdiagrams whose elements satisfy the condition that any 
pairs of them are disjoint or nested. (The empty set is also allowed.) When a forest contains 
the diagram Q itself, it is called full forest. Otherwise it is normal forest. For the calculation 
of g —2 term it is sufficient to consider only normal forests. We denote the set of all possible 
normal forests of a diagram Q by 3(G)- 


B. Forest formula and R-operations 


Assume €5 is a subtraction operator associated with a subdiagram S. The renormalized 

amp,ltu ff of a diagram 5 is obtained frora the u—dizedampMude Me by f ° rest 

formula j 2 j|: 

m s= e < 9i > 

/eSfS) Ue/ 

Here the sum is taken over all forests / of the diagram Q. The order of operation in the 
product is arranged so that the inner subdiagrams are applied first. 

In our approach the subtraction operation is provided as R-operation for performing 
the intermediate renormalization in numerical procedure. Recall that the integrand of the 
subtraction term obtained by R-operation factorizes exactly into the UV-divergent part 
of renormalization constant and the lower-order g — 2 term by construction. This feature 
enables us to apply repeatedly A-operations if there is another UV-divergent subdiagram 
in the forest. 


C. Procedure 

The subtraction term associated with a forest is obtained by the successive A-operations. 
The concrete procedure is given as follows. Suppose the forest / consists of m subdiagrams, 
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/ = {iSi, ..., S m }. They are arranged in such an order that the inner subdiagrams come 
ahead. 

The scaling like Eq. (ESD for a forest / can be defined similarly by introducing the scaling 
parameters {e k } k=l ... m as 

z i — 0( e k) i G <Sk, (92) 

where Sk is the inner-most subdiagram of the forest / which contains the line i. We define 
the UV limit [.9 ({Tj'})]ijv °f a function g{{zj}) of Feynman parameters as the leading term 
of g({zj}) in the successive limits 

kK{%})]uv = lim n ■ ■ ■ iim o({^}). (93) 

6m 61 

Due to the factorization property we can construct the subtraction term corresponding 
to the forest / by the repeated applications of ^-operations. The R-operation for fcth 
subdiagram Sk is applied to the reduced diagram Q /(5i U • • • U5fc_i) which has been obtained 
by shrinking the subdiagrams up to (k — l)th subdiagrams to points. The integrand of the 
subtraction term is obtained in the following way. 


(1) In Eq. 115311 collect all terms which are maximally contracted within the subdiagram S{ 
for i = 1 ,..., 77i. 

(2) Replace U, Bij, Cij and Aj appearing in the integrand by their UV limits, [U](j y , [B r j\ mn 
l C ij\ uv> and [Atluv’ respectively. 

(3) Arrange V in the limit to take the form 


H-H Vs, + ^/(au-us.). 


(94) 


where Sk is a subdiagram obtained from Sk by shrinking all inner subdiagrams to points. 


(4) Attach the overall sign (—l) m . 

The end result of this construction is identical with what was obtained in Refs. QQQ. 
The advantage of the forest approach is that it is readily translatable into computer code, 
which of course is crucial for automation of our entire formulation. 


V. Q-TYPE FEYNMAN DIAGRAMS 


In this section and the next we focus on the particular type of diagrams which have no 
closed lepton loops. We call such a diagram as q-type. The q-type diagram has a simple 
structure, which allows simple identification of various graph-theoretical notions embedded 
in the diagram. The key features relevant for automated scheme of calculations may be 
listed as follows: 

(a) systematic generation of diagrams are easily done. 

(b) a set of independent loops are easily identified. 

(c) subdiagrams relevant for the UV divergence are easily identified. 

They enable us to develop efficient algorithms and implementations. 

In this section the features and © are discussed. A set of algorithms to obtain the 
complete set of topologically independent diagrams is presented. The feature © is related 
to the construction of the topological forms B tJ and U, which provide the building blocks 
of the integrand. The feature (jrj) will be discussed in the next section in relation to the 
subtraction of UV divergences. 

A. Definition and diagram representation 

A q-type self-energy-like diagram of 2 nth order is given by a path V consisting of lepton 
lines emanating from the incoming lepton 0 (y) and terminating at the outgoing lepton ip{x), 
with n photon lines attached to the path at their both ends; there is no closed lepton loop. A 
typical diagram is shown in Fig.0 In QED there is only one type of interaction, the coupling 
of electromagnetic current = - 07^-0 to the gauge potential A^, which is represented in 
Feynman rules as a trivalent vertex, at which a photon line is attached to the lepton line 
path. Here we consider only one-particle irreducible (1PI) diagrams. 

We denote 2 n vertices as Vj (j = 0,..., (2 n — 1)), and adopt the following convention 
throughout the rest of the paper. Every q-type diagram is drawn in such a way that the 
path V passes from the right to the left. The vertices Vj sequentially lie on the path V as 
vq ..., V 2 n-i from the left to the right. A lepton line is denoted as 4 (k — 1,..., (2 n — 1)) 
which runs from a vertex Vk to another vertex Vk-\- A photon line which connects two vertices 
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FIG. 5: An example of q-type diagram. 



Vo Vi v 2 v 3 v 4 v 5 


Vi k and Vj k , is denoted as hk (k — a,b,...), where label k is taken to be an alphabet. The 
photon line is also represented by the pair of two endpoints (v ik , v Jk ). For later convenience, 
a direction of the photon line is chosen by Vi k —> v Jk where the indices are ordered as 4 < jk- 
We also denote the lines collectively as {/}. 

A q-type diagram Q is uniquely specified by the set of n photon lines, i.e., the set of pairs 
of vertices as 

Q = , (95) 

where 4, j i, ... ,i n , j n take values in {0,..., (2 n — 1)} exclusively. To avoid the ambiguity of 
representation, we further impose the following condition: 


4 < jk , for each pair (v ik ,v jk ), 

4 < 4 q between two pairs k < k'. 


(96) 


B. Circuits and loop matrix 


Consider the chain diagram associated with a q-type diagram. A circuit of the graph is 
a self-nonintersecting closed path along the graph. The maximal independent set of them 
defines a fundamental set of circuits, which provides a complete basis of closed loops of the 
graph. The fundamental set of circuits of a graph is crucial in the calculation of U and B tJ . 

In the case of q-type diagram, a path C r composed of a photon line h r = (v ir , v Jr ) and 
lepton lines U r+ 1 , ..., l Jr that lie between two endpoints of the photon line form a circuit. 
The direction of the circuit can be taken naturally along those of lepton and photon lines. 
The fundamental set is thus chosen by the set of paths {C r } r=1 n . 

The loop matrix £j t c r is given almost trivially. Once ^ h c r is known, U and are calcu¬ 
lated according to the formula (HD, (HD, and m given in Section El or (BSD, (ED, and 
(BgP in Appendix E 
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C. Time-reversal symmetry 


By using the time-reversal symmetry of QED, we can further relate distinct diagrams 
to each other, to reduce independent set of diagrams. Two q-type Feynman diagrams are 
equivalent in time-reversal and give rise to the same contribution to anomalous magnetic 
moment if and only if they are the images of each other under the reversal of the directions 
of all lepton lines in V. 

For a time-reversal equivalent pair of diagrams, it is sufficient to evaluate either of the two. 
The asymmetric diagrams in time-reversal dominate at higher orders, which implies that the 
number of distinct diagrams to be evaluated is cut down to almost half by considering time- 
reversal symmetry. 

D. Algorithms 

Suppose 2n vertices are placed on a lepton line path V. Then the complete set of topo¬ 
logically distinct q-type diagrams of 2ntli order is obtained as follows: 

Step 1. Connect a pair of vertices by a photon line in every possible way. 

Step 2. Pick only 1PI diagrams and discard others. 

Step 3. Drop either of the pair of equivalent diagrams under time-reversal. 

Step 1 is a process to list up all possible ways to make n pairs out of the vertices 
{0,..., (2 n — 1)}. A procedure to make k pairs from 2k elements is given as follows; the 
elements are assumed to be ordered in a line. 

Pick an element at the left end of the line, and another one from the rest to form 
a pair. Repeat the process to the remaining 2k —2 elements to make k — 1 pairs 
until there is no element left. 

By considering 2k — 1 ways to make a pair, we can generate all possible pairings recursively. 
The total number of ways is {2k — 1)!!. 

Next we go to step 2. The diagram corresponding to a pairing generated above may not 
satisfy the 1PI condition. Since a q-type diagram is connected by lepton lines, it is sufficient 
to check whether it stays connected when one of the lepton lines is eliminated. 
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A q-type diagram is 1PI if and only if for each lepton line Ik, there exists at least 
one photon line that steps over the lepton line, i.e ., two end points of photon 
line, (vi,Vj), satisfy ry < Vk-i and Vj > Vk simultaneously. 

The diagrams that do not match the condition shall be discarded. 

The time-reversal operation in step 3 is done by substituting the index k of vertex by 
2n —1~ k. A q-type diagram Q is mapped to Q' by the substitution of indices followed by the 
reshuffling of pairs to satisfy the conventions ©• If Q is invariant under the time-reversal, 
it should be kept with the symmetry factor one. Otherwise, either of Q or Q' should be kept 
with the symmetry factor two; we adopt the rule that the diagram Q is chosen when the 
lexicographical order of the patterns of indices representing the diagrams Q is ahead of Q'. 

E. Number of diagrams 

Based on the above consideration, the number of 1PI q-type diagrams J\f n of n loops is 
given recursively by the following relation (disregarding time-reversal symmetry): 

U n = (2n-l)\\- n^. ( 97 ) 

where denotes the set of ordered partitions of n (i.e. (1,2) and (2,1) should be distin¬ 
guished) . 

Table ID shows the number N n of independent q-type diagrams for n < 7 as well as that of 
symmetric ones and that of asymmetric ones under time-reversal. (J\f n = N n sym + 2V n asym .) 
It also demonstrates that the incorporation of time-reversal symmetry efficiently reduces the 
number of independent diagrams to be evaluated (N n sym + A n asym ) at higher orders. 

VI. UV DIVERGENCE STRUCTURE OF Q-TYPE DIAGRAMS 

For a given Feynman diagram, it is required in the renormalization process to pick up 
all the 1PI subdiagrams that have overall ultraviolet (UV) divergence. They are referred 
to as UV divergent 1PI subdiagrams. By the power counting of superficial divergence, 
there are two types of UV divergent subdiagrams in q-type diagrams in QED, namely, the 
lepton self-energy-like subdiagram, and the vertex subdiagram. Every subtraction term 
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TABLE I: The number of independent 1PI diagrams of q-type, N n . with n loops. 


n 

N n 

N sym 

iy n 

at asym 
ly n 

1 

1 

l 

0 

2 

2 

2 

0 

3 

8 

6 

2 

4 

47 

20 

27 

5 

389 

72 

317 

6 

4226 

290 

3936 

7 

55804 

1198 

54606 


in the subtractive renormalization procedure corresponds to the Zimmermann’s forest, a 
combination of subdiagrams whose loop momenta go to infinity. 

For a q-type diagrams, the above prescription is implemented quit simply, reflecting the 
graph theoretical properties of the diagram. In this section we discuss the UV structure of 
the q-type diagrams and describe an algorithm to compose subdiagrams and forests of the 
diagram. 

A. UV divergent subdiagrams 

A subdiagram relevant for the UV divergence is either of the self-energy type or of the 
vertex type. For a q-type diagram in which all vertices are located on the lepton line path 
V a subset of vertices are denoted by one or more segments of the path. Thus a subdiagram 
of these types of a q-type diagram corresponds to a single segment of the path, and it is 
specified by the indices of two end-point vertices. 

Therefore, to obtain all the divergent subdiagrams of a q-type diagram we have only to 
find every possible pair of indices [i,j\, 0 < i < j < (2 n — 1) that satisfies the following two 
conditions: 

The subdiagram corresponding to the segment [i,j] is classified into the self¬ 
energy-like type or the vertex type. The number of ‘floating’ photon line (only 
either one of the two endpoints of the photon line lies on the segment [i,j]) is 
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zero (for self-energy-like subdiagram) or one (for vertex subdiagram). 

And, 

It is one-particle irreducible, i.e ., it stays connected when any one of the lepton 
lines that belong to the subdiagram is eliminated. 

The second condition is satisfied when for every lepton line, U+ 1, ■ ■ ■ ,lj, there is at least one 
photon line that belongs to the subdiagram (both endpoints of it lie between w* and Vj) 
which steps over the lepton line. The photon line hk = (vj k , Vj k ) steps owe?’the lepton line 
l s = (w s _i, w s ) when i < i k < s — 1 and s < jk < j simultaneously. 

It is noted that subdiagrams of q-type diagrams also do not contain lepton loops. The 
residual diagram of a q-type diagram which is obtained by shrinking the subdiagram to 
a point is again of q-type. Recall that it is related to lower-order g — 2 term. The UV 
subtraction procedure is closed within the q-type diagrams. 

B. Forests 

To begin with, we define the inclusion relation between two subdiagrams, S a and Sb- 

disjoint if S a and Sb do not share any vertices nor lines, i.e. S a f] Sb = 0. 

overlapping if S a and Sb share some vertices and lines though one is not completely in¬ 
cluded in the other. 

nested if S a (or Sb) is a subset of the other, i.e. S a C Sb or S a D Sb- 

For the q-type diagrams, the inclusion relation is mapped to that of two segments. The 
relation between two subdiagrams represented by S a = [i a ,j a ] and Sb = [ib,jb] is one of the 
following (assuming that i a < %)■ 

disjoint if j a < ib, 
overlapping if i a <i b < j a < jb, 


nested if i a < ib and jb < j a - 


1a Ja 1b Jb 

1 , \ Ja __ 


lb 


Jb 


1a ^b Jb Ja 
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FIG. 6: A forest composed of nested subdiagrams 71,... ,74 (left), and the corresponding cascade 
structure (right). 



72 74 



A forest is defined as such a set of subdiagrams that any two of its elements are not 
overlapping with each other; they are disjoint or nested. Since we are currently interested 
in the magnetic form factor, it is sufficient to consider only the ‘normal’ forests which do 
not contain the diagram Q itself. On the other hand, a forest that contains Q is called ‘full’ 
forest. 

A complete set of forests of a diagram is generated by finding all the combinations of the 
subdiagrams, and discarding those which contain the overlapping subdiagrams. (This partic¬ 
ular procedure is not restricted to the q-type diagrams.) A cascade structure of subdiagrams 
of a forest is reproduced by referring to the inclusion relation between subdiagrams. This 
information is required during the subtraction process which is performed for the divergence 
corresponding to the inner subdiagrams first. 

VII. AUTOMATED FLOW OF CALCULATION 

In this section we present the flow of the process to generate the numerical integration 
code for evaluating an individual diagram from its representation indicated by the rectangu¬ 
lar box at the upper-left corner of Fig. 0 The provided information enables us to construct 
the amplitude in the form of Feynman parametric integrals in terms of building blocks, U, 
Bij, scalar currents, and so forth. This follows exactly the pattern developed for the sixth- 
and eighth-order cases [hi fici ] . Next, the ultraviolet divergence is treated via A-operation 
which identifies and subtracts the most divergent part of the original integral, corresponding 
to a specific UV limit. The treatment of the whole divergence structure is organized by the 
Zimmermann’s forest formula. The infrared divergence remaining in the individual diagram 
should also be subtracted away, though this article does not cover this subject. Finally, the 
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FIG. 7: Flow of process to generate the numerical integration codes from the diagram representa¬ 
tion. 



(intermediate-) renormalized amplitude constructed from the original amplitude and the 
set of subtraction terms is turned into a FORTRAN code, which is readily processed by 
the numerical integration system such as VEGAS 1-25] , an adaptive Monte-Carlo integration 
routine. 

A. Diagram generation 

We begin by generating a complete set of topologically distinct q-type diagrams of a given 
order according to the algorithm in Section IV D1 The implementation of the algorithm is 
achieved in C++. Each diagram is expressed by a single-line representation which describes 
the pattern of pairings of vertices by the photon propagators. The diagrams are then named 
after a certain convention and stored in a plain text file. All subsequent steps refer to this 
hie for the diagram data. 


36 











































































Names and forms of all relevant diagrams of the sixth- and the eighth-orders are listed 
in Refs. [7, 12 I |. We adopt the following convention for the tenth-order diagrams. The W-T 
summed diagrams are classified into two groups, one of which is time-reversal symmetric (72 
diagrams) and the other is asymmetric (317 diagrams). They are sorted in a lexicographical 
order within each group, and then given serial numbers with a prefix “X” which stands for 
the tenth-order in a Roman numeral, first for the symmetric ones (X001, ..., X072), then 
for the asymmetric ones (X073, ..., X389). 


B. Subdiagram search and forest construction 

All UV divergent subdiagrams of a q-type diagram are identified according to the al¬ 
gorithms in Section IVIA1 Then all forests are constructed according to the description in 
Section IViBI bv generating every possible combination of subdiagrams which does not con¬ 
tain any overlapping pairs. The inclusion relation between subdiagrams is examined prior 
to this step. 

The cascade structure among the subdiagrams of a forest is also identified and stored in a 
tree form. It is, however, not mandatory since the order of successive subtraction operations 
of a forest is automatically respected if we adopt the regulation that the diagrams of smaller 
sizes are processed first. It is because the sizes of subdiagrams 5) and Sj satisfy S', < Sj 
whenever S t is contained in Sj. 

The implementation is carried out in both Perl and C++. To demonstrate the fast 
algorithms and implementations we generated the diagrams and their forests up to 14th 
order, which took less than 10 minutes on an ordinary PC. 

C. Constructing unrenormalized integrand 


A single-line representation of a q-type diagram may be directly translated into the form 
of unrenormalized integrand. Recall that the integrand F of a q-type diagram is given by 
Eq. (jSJ): 


F = (]pi + mi) • • • 7 /i2n_1 (^ 2 n-i + m 2n - 1 ) 7^ 2n x JJ 9m 


k^ik ’ 


(98) 


k=1 
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where the diagram-specific product n g, Mk lljk in this case is determined by the pairing pattern 


{ («ii 5 Vji ) i {Pi 2 , Vj 2 ))■•■} • 


(99) 


The basic form of integrand is common to all q-type diagrams, where we have only to make 
permutation of indices to construct the integrand of a particular diagram. Therefore we can 
make full use of a template to perform this step. 

We implemented this step as a Perl program, which translates the single-line representa¬ 
tion of the diagram into an explicit form of F and put it into a template of FORM program, 
ft is then processed by FORM to perform the analytic integration over all loop momenta, 
the trace calculation, contractions of D operators, and other algebraic manipulations, which 
yields the unrenormalized integrand expressed as a polynomial of the building blocks, U, B^, 
Cij and Aj, integrated over Feynman parameters. This follows exactly the method developed 
for the sixth- and eighth-order cases flil lbf| . The output is in the FORM-readable form 
for the subsequent steps of deriving various UV limits, as well as in the form of FORTRAN 
code. 

D. Constructing building blocks 


The building blocks of infegcand, U and B,„ ace defined fan. the undec.y.ng topo- 
logical structure of the diagram called chain diagram 11 -1 121| . First we identify chains 
and chain variables z a = JV ga A- The fundamental set of circuits {C r }r=i....,n of the dia¬ 
gram is identified according to the specification in Section IV Bl and the loop matrix £ a> c r 
is constructed accordingly. Once £ a ,c r is known the building blocks U, are obtained as 
homogeneous polynomials of { z Q } by Eqs. USED, (El, and (ED- Other building blocks, C tJ 
(C^ = Cij/U ), Aj, and V are also constructed in terms of U, B^ and Feynman parameters 
{zi} by Eqs. (jHIlj) . ED, ED, an d ED : 


1 


4 ? — JJ Zi ^ij ’ 

i 

\ = Zi — G , G = Zi Ai , 


( 100 ) 

( 101 ) 


where the sum runs over all lepton lines. 

The calculation of U, Bij and C tJ involves algebraic manipulations such as determinants 
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and cofactors of matrices whose elements are polynomials of Feynman parameters {z a }. We 
implemented this step in three ways. 

i) The algebraic manipulations are performed by MAPLE. We identify the loop matrix 
from the diagram representation and prepare a MAPLE program to calculate U, B^, 
and Cij according to Eqs. (EH, m, ED, and (EH). The output is in the FORM- 
readable form for the subsequent operations. The FORTRAN code is also generated 
from it via FORM. 

ii) We developed concise formulae, (IB 51) and (1B6D . to provide the coefficients of the U- 
function directly from the loop matrix. The coefficients u PlP2 ,„ Prn are defined by 

v = Y, »*.* < < ■ ■ ■ <” . ( 102 ) 

{Pi} 

m 

for every possible combination of {pi}, where pi takes 0 or 1 and pi = n. m is 

i=1 

the total number of chains. Similarly for Bij and Cij by the formulae described in 
Appendix El They are implemented in C++. 

iii) The algebraic manipulations are also handled in C++ by constructing proper data 
structures (or classes ). We developed a simple polynomial class and implemented the 
calculation of U and B l3 according to Eqs. El, El, and ED- We also developed 
another version with the help of GiNaC ^8|, an algebraic manipulation library in C++. 

The current version of automation system relies on the implementation 13) for no particular 
reason. The scalar currents Aj and the function V are constructed from U and B %3 according 
to Eqs. (11001) and CHD- 

E. Constructing UV subtraction terms 

The UV divergences of a diagram are identified as Zimmermann’s forests which are 
the combinations of UV divergent subdiagrams. We construct a set of UV subtraction 
terms, each of which corresponds to a particular forest, by the successive applications of 
R-operations to the unrenormalized integrand. 

The subtraction term is constructed by the following three steps. 
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1. Find the UV limit of the building blocks. 

2. Find the UV limit of the integrand. 

3. Modify the UV limit of V-function in the denominator to satisfy the fac- 
torizability requirement of subtraction terms. 

Those steps are achieved by simple power counting applied to the original (unrenormalized) 
integrand and building blocks, without referring to any lower-order constructs. The whole 
implementation is done in Perl with the help of MAPLE and FORM for symbolic manip¬ 
ulations. This follows exactly the scheme developed for the sixth- and eighth-order cases 

m. 


1. UV limit of building blocks 

The explicit form of building blocks, U, Bij, and C l3 , in the UV limit (EHD related to a 
subdiagram S is given as the leading term in power series expansion by e under the rescaling 
of the Feynman parameters as 

Zi —* e Zi , % E S . (103) 

The procedure is implemented as a MAPLE or FORM program, in which the scaling rules 
are generated from the information of the subdiagram. 

For a forest consisting of more than one subdiagrams the above procedure is successively 
applied with each subdiagram The order of operations is determined referring to the 
cascade structure of the forest so that the inner subdiagrams are applied first. 

The UV limit of scalar current Aj is constructed from the UV limits of U and B i3 . 


2. UV limit of the integrand 

According to the formulation in Section IIII A1 the UV divergent part of the integrand in 
the UV limit (15511 is derived from the most contracted terms within the subdiagram S. This 
part is simply extracted by counting the number of B i3 with i,j E S in each term of the 
unrenormalized integrand. 

The procedure is implemented as a FORM program, which reads the expression of inte¬ 
grand constructed in Section IV1I Cl and picks up the terms which are the products of the 
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specified number of Bij whose indices belong to the subdiagram S. For a forest consisting 
of more than one subdiagrams, the above counting are applied successively according to the 
order with the inner subdiagrams first. The FORM program is generated referring to the 
forest data. 

3. UV limit of V-function in the denominator 

The UV limit of V-function in the denominator Vg/s is replaced as follows: 

Vg/s —> Vs + Vg/s (104) 

in the step (3) of if-operation to guarantee the factorization property. At a glance this 
operation might require the explicit construction of V-functions of lower order diagrams, 
Vs and Vg/s individually. However, it turns out that this can be achieved by adopting the 
following rule: in the construction of V by Eq. dUD the scalar currents Aj with j G S 
should be replaced by [A,-] 5 , where [A,-] 5 is given by 

(a) dropping all terms containing B %] with i G S and j G G/S, 

(b) replacing other B tJ and U by their UV limits. 

Thus the replacement of V-function is also accomplished solely by the limit operations from 
the original building blocks. 

F. Symbolic expressions of subtraction terms 

The subtraction term has a symbolic expression in terms of the product of renormalization 
constants and lower order magnetic moment part, each term of which is related to the 
particular structure of the corresponding forest. The identification of the symbolic expression 
is achieved by pattern matching based on the rule set. 

We prepare the rule set for recognizing the particular pattern of subdiagrams (after 
shrinking the inner subdiagrams to points) and identify the form of the expressions. The 
whole implementation is done in Perl and the rule set is also generated automatically from 
the basic representation of the self-energy-likc diagrams. 

The symbolic expression is better suited for human recognition. It will also be relevant 
when we perform the residual renormalization. 
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G. Controlling the whole steps 


Each step of code generation is achieved by separate Perl programs with the help of 
MAPLE and FORM, while the flow of the whole process is governed by a shell script. It 
takes the name of the diagram as an input and performs the following operations: 

a) Pick out the corresponding expression of the diagram from data file. 

b) Construct each component of the integration code. 

c) Gather up the FORTRAN codes in the end. 

The whole process of code generation for each W-T summed diagram of tenth order takes 
10-20 minutes on an ordinary PC. 


VIII. CONCLUSION AND DISCUSSIONS 

In this paper we presented an automated scheme of code generation for evaluating higher 
order QED corrections of the lepton anomalous magnetic moment by means of numerical 
method. We constructed an algorithm and concrete procedure to obtain UV-hnite ampli¬ 
tudes for a particular set of diagrams without lepton loops, which we call q-type diagrams. 
Our current concern is the tenth-order corrections, though, the scheme itself is applicable 
to an arbitrary order. 

We implemented our procedure as a set of Perl programs with the help of symbolic 
manipulation systems, FORM and MAPLE. From a single-line representation of a diagram 
it generates numerical integration codes in FORTRAN, which are ready to be processed by 
VEGAS, an adaptive Monte-Carlo integration routine. 

The programs have been tested for lower-order diagrams and confirmed that they re¬ 
produce the codes for the sixth order and eighth order diagrams previously constructed. 
They are now being applied to tenth-order diagrams. At present, the diagrams which have 
only vertex renormalization were processed and test runs were performed. Those diagrams 
corresponds to 2232 vertex diagrams among 6354 q-type diagrams of tenth-order. Crude 
evaluation showed no sign of divergent behavior, which confirms that our scheme is working 
well. They are currently put to production runs. 
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The remaining 4122 diagrams have not only UV divergent self-energy subdiagrams but 
also infrared (IR) divergences. The simplest way to deal with the IR problem is to give a 
small mass A to photons, which requires no further work on the automating code, and is 
being pursued as the first step. 

This scheme has thus far been tested successfully with the sixth-order q-type diagrams, 
and reproduced the analytic result after proper treatment of residual renormalization terms. 
For the eighth-order q-type diagrams it seems successful so far except for a few diagrams 
which suffer from more severe IR divergences than logarithmic. One remedy is to subtract 
full part of renormalization term by taking properly into account the effect of self-mass 
term of the leptons. This modification will be implemented within a slight extension to the 
current automation code and is being worked out. 

To obtain a result independent of A it is necessary to incorporate IR subtraction terms in a 
manner similar to that of UV counterterms. The subsequent step of residual renormalization 
is our next issue. 

Our scheme has been elucidated for the q-type diagrams in this paper, though the formal¬ 
ism is not restricted to that type of diagrams. The practical algorithm for the construction 
of building blocks, that are related to the underlying topological structure of the diagram, 
and the identification of UV divergent subdiagrams rely on the particular properties of the 
q-types. However, they can be extended to incorporate more general cases. We will then 
have a fully automatic scheme for evaluating QED diagrams of lepton g — 2. 
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The numerical calculation has been performed on the RIKEN Super Combined Cluster 
System (RSCC). 

APPENDIX A: CLASSIFICATION OF DIAGRAMS CONTRIBUTING TO Aj 10) 

There are 12672 vertex-type Feynman diagrams at the tenth order. We classify them into 
six sets according to the type of virtual lepton loop(s) and how they appear in a Feynman 
diagram. Every figure in this appendix should be supplied with one external vertex in all 
non-trivial places of the internal lepton lines. 

All diagrams in the sets I, II, III and IV are obtained by inserting vacuum-polarization 
and/or light-by-light-scattering subdiagrams of appropriate orders into lower-order q-type 
diagrams. 












FIG. 8: Set I. There are 208 Feynman diagrams in this set. 

Set I consists of 208 diagrams of the form shown in Fig. |H1 Each of these diagrams is 
obtained by inserting vacuum-polarization or light-by-light-scattering subdiagrams into a 
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q-type diagram of the second order. Set I can be classified into ten gauge-invariant subsets. 


1(a) A diagram contains four vacuum-polarizations of the second order. 

1(b) Each diagram contains a fourth-order vacuum-polarization and two vacuum- 
polarizations of the second order. 

1(c) Each diagram contains two vacuum-polarizations of the fourth order. 

1(d) Each diagram contains a second-order vacuum-polarization and a sixth-order vacuum- 
polarization which consists of two lepton loops. 

1(e) Each diagram contains a second-order vacuum-polarization and a sixth-order vacuum- 
polarization which consists of a single lepton loop. 

1(f) Each diagram contains an eighth-order vacuum-polarization which consists of three 
lepton loops. 

!(g) Each diagram contains an eighth-order vacuum-polarization which contains a fourth- 
order vacuum-polarization as its subdiagram. 

1(h) Each diagram contains an eighth-order vacuum-polarization which contains a second- 
order vacuum polarization as its subdiagram. 

I(i) Each diagram contains an eighth-order vacuum-polarization which consists of a single 
lepton loop. 

I(j) Each diagram contains an eighth-order vacuum-polarization which consists of two light- 
by-light-scattering subdiagrams. 

Set If consists of 600 diagrams of the form shown in Fig. El each of which is obtained 

by inserting vacuum-polarization subdiagrams and/or a light-by-light-scattering subdiagram 

into a q-type diagram of the fourth order. Set 11 can be further classified into six gauge- 

invariant subsets. 

11 (a) Each diagram contains three vacuum-polarizations of the second order. 

11(b) Each diagram contains a second-order vacuum-polarization and a fourth-order 
vacuum-polarization. 
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FIG. 9: Set II. There are 600 Feynman diagrams in this set. 

11(c) Each diagram contains a sixth-order vacuum-polarization which contains an internal 
lepton loop. 

11(d) Each diagram contains a sixth-order vacuum-polarization which consists of only one 
lepton loop. 

11 (e) Each diagram contains a sixth-order light-by-light-scattering subdiagram. 

11(f) Each diagram contains a fourth-order light-by-light-scattering subdiagram and a 
second-order vacuum-polarization. 





FIG. 10: Set III. There are 1140 diagrams in this set. 

Set III consists of 1140 diagrams of the form shown in Fig.^B each of which is obtained by 
inserting a vacuum-polarization subdiagram and/or a light-by-light-scattering subdiagram 
into a sixth-order q-type diagram. Set III can be further classified into three gauge invariant 
subsets. 
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ni(a) Each diagram contains two vacnnm-polarizations of the second order. 

111(b) Each diagram contains a fourth-order vacuum-polarization. 

ni(c) Each diagram contains a fourth-order light-by-light-scattering subdiagram. 



FIG. 11: Set IV. There are 2072 Feynman diagrams in this set. 


Set IV consists of 2072 diagrams of the form shown in Fig. each of which is obtained 
by inserting a second-order vacuum-polarization into an eighth-order q-type diagram. 

Set V in Fig. IT?1 consists of 6354 vertex diagrams of the q-type of the tenth order. 



FIG. 12: Set V. There are 6354 Feynman diagrams in this set. 

Set VI consists of 2298 Feynman diagrams of the form shown in Fig. El They contain a 
light-by-light-scattering subdiagram, one of whose photon lines is supposed to be external. 
We call them as l-by-l-type diagrams hereafter. This set is further classified into eleven gauge 
invariant subsets. Each subset of diagrams also includes radiative corrections of respective 
types except for the subset VI(k). 

VI(a) Each diagram is obtained by inserting two vacuum-polarizations of the second order 
into a 1-by-l-type diagram. 

VI(b) Each diagram is obtained by inserting a fourth-order vacuum-polarization into a 
1 -by-l-type diagram. 

VI(c) Each diagram is obtained by inserting a second-order vacuum-polarization into one 
of virtual photon lines coming out of the lepton loop of the 1-by-l-type diagram and 
attaching one internal photon line to the open lepton line. 
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FIG. 13: Set VI. There are 2298 Feynman diagrams in this set. 

VI(d) Each diagram is obtained by attaching two internal photon lines to the open lepton 
line of a 1-by-l-type diagram. 

VI(e) Each diagram is obtained by attaching an internal photon line with a second-order 
vacuum-polarization inserted to the open lepton line of a 1-by-l-type diagram. 

vi(f) Each diagram is obtained by first attaching an internal photon line within the lepton 
loop of a 1-by-l-type diagram and then inserting a second-order vacuum-polarization 
into one of the three photon lines which connect the lepton loop and the open lepton 
line. 
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VI(g) Each diagram is obtained by first attaching an internal photon line within the lepton 
loop of a 1-by-l-type diagram and then attaching an internal photon line to the open 
lepton line. 


VI(h) Each diagram is obtained by attaching two internal photon lines within the lepton 


loop of a 1-by-l-type diagram. 

VI(i) Each diagram is obtained by attaching an internal photon line with a second-order 
vacuum-polarization inserted to the lepton loop of a 1-by-l-type diagram. 

via) Each diagram is obtained by inserting a light-by-light-scattering subdiagram of the 
fourth order into a 1-by-l-type diagram 

VI(k) Each diagram contains a light-by-light-scattering amplitude with six photon legs. 

APPENDIX B: ALTERNATE FORMULAE FOR U AND B POLYNOMIALS 

Definitions (BUI) and (Uffl) are useful for the programs like MATHEMATICA and MAPLE. 
However, they turned out to be clumsy for developing programs by languages such as C++ 
and FORTRAN. This is why it is useful to find alternate formulae for U and B polynomials. 
Our strategy is to directly provide compact expression for each term which is homogeneous 
polynomial of Feynman parameters {zi}. 

1. Concise formula for U 

Let us first consider the U polynomial. We introduce the circuit matrix characterized by 
the chain indices a instead of the line indices j by 


£,a,r £,j, r (j £ Oi) . 


(Bl) 


This is unambiguously defined in our convention for orientations of chains and lines. Then, 
Eq. (1291) can be rewritten in the form which expresses that U and U rs are determined solely 
by the structure of chains: 



(B2) 


all a 
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where w a is defined by w a 


J2ica z i- Now > U of Eq. m can be expanded as 


U = det U rs 

1 <r, s<n 

= ^ ^ w ai £<* 1 , 1 ' ' ' 'y ^ w oin £>an,n 
ai On 

x ^ ] ^( cr ) Coi, O-(I) ■ ■ ■ £ct n , a(n) 1 (B3) 

rrG&n 

where & n denotes the permutation group of degree n and e(a) = ±1 is the signature of 
cr G & n . The right-hand side of the last equality shows that the terms with a r = a s for 
r ^ s vanishes. Thus U is a homogeneous polynomial of degree n, where each monomial 
can be at most linear with respect to each w a . Such a monomial is characterized by a 
combination {cti, • • • , a n } whose elements are picked up from the sets {1, • • • , 3n — 3}. 
They can be ordered as ay < ■ ■ • < a n taking the permutation of the indices attached to the 
circuits in Eq. (IB3D into account. Then, Eq. (IB3I) becomes 

U = w ai ---w an 

l<ai<---<an<3n—3 

* ^ 1 £ai,<T'(l) ’ ’ ’ Ca„,a’(n) 

a'G&n 

X ^ ] ^( cr ) £ai,a(a'(l)) ’ ’ ’ £a n ,o-(<r'(n)) • (B4) 

c tS6„ 

We replace the sum over a € & n by the sum over cr" = a o a 1 G & n . Then, noting 
e(cr") = £(a / )e(cr), the quantities appearing in the second line and the third line on the 
right-hand side of Eq. (1B4I) turn out to factorize and coincide with each other. Thus, we get 
the following formula for U: 


U= ^2 w ai -■-w an (A(a 1 , ■ ■ ■ , a n )) 2 , (B5) 

l<ai<---<an<3n—3 

with 


4(«i, 


(y,r 


) — ^ 1 £ai,<r(l) ■ ■ ■ £a 

cr£&n 


r{n) 


(B6) 


2. Concise formula for B a p 

We also find a convenient formula for B polynomials following a similar manipulation as 
above. We start with the expression 

n 

B a0 = uY J Ur{U~ l ) rs ^s, (B7) 

r, s =1 


50 


showing that they are also determined solely by the associated chain diagram. By deferring 
its derivation below, we obtain a compact formula 






><A(a, i &n— l) A(^/3 ^ ^1? ? l) 5 


(B8) 


with A(a, a±, • • • , a n _i) given in Eq. (IB6I) . For a = f3, there is a particular equation which 
follows from Eqs. dm) and m 0 , 

dU 


B n rv 


DWn 


Previously most hand calculation has been done using the Nakanishi formula 

Ba 0 = 'y ^ Ca, c C/3, c Ug/ C , 


(B9) 


(BIO) 


where the sum runs over all possible circuits (not necessarily limited to those of the funda¬ 
mental set of circuits), Ca,c = (1> — 1 ; 0) is a projector to whether the chain a runs (along, 
against, outside of) the loop c, and Ug/ C is the U polynomial of the reduced diagram Q /c 
obtained by shrinking the circuit c and the vertices on it to a single vertex. The formula 
(IB 101) requires a topological manipulation to pick up all circuits in a graph, while the formula 
enables us to calculate B polynomials algebraically. 


Derivation of Eq. (IB81) . 

We derive the formula dBSl) for B a p. For that purpose, we introduce a set of functions 
{/ r } r=1 ... n on the domain {1, • • • , n — 1} by 


fr(x) = 


X 


X + 1 


(1 < x < r — 1) 
(r < x < n — 1) 


(Bll) 


Then, the elements of the (n — 1) x (n — 1) minor matrix Uf s obtained by eliminating the 
stli column and the rth raw from {U sr } s r=1 n are given by (also using the symmetric 
property of U rs ) 

(Ufs) xy = Uf r ( X ),f s (y) ■ (B12) 

Inserting Eq. (IB 1 2|) into 


u (U' 1 ),, = (- 1 ) r+s „ d s „((W,W, 

l<a;, 2 /<(n—1) 


(B13) 


51 








the part U(U 1 ) rs in Eq. (IB7D can be written 


as 


U(U 




< 1) + £ ( a )Uf r (l),f a (<j(l))---Ufr( n -i)J s (cr(n- 1)) 

cr £& n -1 

Hr +S £ W Ct 1 £ai,fr( 1) E ®On-l ^d„-l,/ r (fl-l) 

«i a n _i 

X y ' £ ( a ) £cn,f s (a(l)) ' ' ' £a n _i,/ s (o-(n-l)) 

cre©n-l 

^ — 1 ) ^ 1 ' ‘ ' w a n _i £«i, f r (l) ' ' ' £a n -i, f r (n-l) 

a r'¥ za s' 

X y 1 £ ( C7 ) ^Ol,/s(o-(l)) ■ ' ' ^CKn- 1 ,/s(cr(n-l)) 

0-e©n-l 

i-ir- v w^-w^ 

Qi<"-<a„-i 

X ^ q < 7 '( 1 )>/'-( 1 ) ■' '^ a a'(n-l)> fr{n- 1) 

O-'eSn-l 

X £ ( <J ) ^'(p./sM 1 )) ‘ ‘ ' ^ a a'(n-l)Js(cr{n-l)) 

(7$. 

i-i r- E 

<■■■«*., 

y 1 £ai, ■ ■ ■ £a n -l,fr(<r'(n—l)) 

O-'GSn-l 

X ^ 1 £ai,/ s ((<7ocr')(l)) ‘ ’ ’ £a n -i, f s ((cro<7')(n-l)) 


<tG& „ 




cre&n-l 


(-i) ,+ ‘ E 


ai<-<a„_i 


X I £ ( a ) Cai,/r(o-(l)) ’ ‘ ‘ £an-l,/r(o-(n-l)) 

y( 7 GS„_l 

X I £ai,/ s (<r(l)) ■ ■ ■ ^a n _i,/ s ((T(n-l)) 

\(jGS„_i 


By inserting this expression into Eq. m , we get 




where 


(B14) 


B a fj ^ ^ H'n i ■ ■ m w anl B(a, ol i, j Qfn— i) B(y(3^ Oi\, j oi n — i), (B15) 


-B(a;; Op, ••• , Cpi-l) — ^ 1) £a,r ^ ^ e (°") £ai,/ r ((x(l)) ■ • • £an-l, /r (<*(«—1)) • (B16) 


r=l 


CTGSn-l 
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The remained task is to demonstrate that this B(a\ cq, • • • , ct n _i) coincides with 
A(a, a± ■ ■ ■ , a n - 1 ) appearing in the coefficient (IB6I) of each monomial of U. For that pur¬ 
pose, the summation over all permutations of the circuits in a fundamental set is replaced 
by that over all permutations of indices distinguishing the chains; 

n 

B(° «1, ) ^n-l) = £ ( ff ) ^(i),/r(l) ' "^(„-i),/r(n-l) 

r= 1 cr£6 n -1 

n 

= ( — 1 ) £ ( cr ) £a CT (i), 1 ' ' ' Ca CT ( T ._i),r-l £a,r 

r =1 cre© n _i 

^^•cr(r)5^Vl 1), n 

— ^ 1 £ ( CT ) £a,<r(l) Cai,cr(2) ' ’ ' £a n _i,<x(n) 

<r£&n 

= A(a, ai,---, a n - 1). (B17) 


In the above, we use the fact that a sequence of permutations 

[c^cr(l)) j OL a (r— 1); Cp Chr(r) ) j C^cr(n—1)] 

I * [«, Cfc(r(l); 5 ^cr(r — 1) ? ^cr(r) i ? ^a(n— 1)] 
f — > [cp Oil, • • • , « n _i] , 


(B18) 


for a e ©n-i, combines to form all possible permutations of degree n and each step produces 
a signature (—l) r+1 and e(cr) respectively. Therefore, we obtain the desired result (IBSI) . 


3. Concise formula for (7, 




fa — j T Cji j 2 1 


1 

U 

1 

u. 


A (7-polynomial is a function characterized by a pair (j i, j 2 ) of two indices of lepton lines 

(B19) 

Cjij 2 — jf ^ ^fci z k 2 {B kl n Bk 2 j 2 ~ B kl n j (B20) 

k\<k 2 

where the summation ranges over all pairs (k\, k 2 ) of the indices of lepton lines. The reason 
why we call this function as a polynomial will be clarified shortly. From its definition, 
Cj x j 2 = ~Cj 2 j 1 . Thus, we can assume that j\ < j 2 without loss of generality. 

The expression of Cj 1 j 2 needs a little bit care. The factor — prior to the summation 

might imply that C n n develops a singularity like — when U —> 0. If it were the case, the 
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maximally contracted terms containing Cj lJ2 would dominate over the other terms in the 
UV-limit. This should not be the case since the degree of ultraviolet singularities of the terms 
containing C'-polynomials, would be higher than that expected from the original expression 
of Feynman amplitude. A closer examination shows that one factor of U factorizes out from 
[B' kl jl B' k2j2 — B' ki j2 B' k2ji ) in Eq. (1B20D . Thus, C jlj2 are actually polynomials, and the 
maximally contracted terms of C'-polynomials have the same degree of singularity as that 
of the other maximally contracted terms. The factorization of U could be realized at the 
numerical level with use of Eq. (IB201) . However, to avoid round-off errors, it is desirable to 
obtain an analytic expression that calculates C n32 as a polynomial. Below, we write down 
such an expression that also enables us to control the UV limit of the C'-polynomials. 

Since B' k ■ differs from B k j when j — k (See Eq. (ffiflB . we define the quantity 

c( h , fc 2 ; ji, :h) = jj z kl z k2 (B ki n B' k2 j2 - B' kl j2 B' k2 jl ) . (B21) 

and distinguish the six cases; 

I: ji = h < Jj‘2 = k 2 

c {jli 32, jl, 32 ) = Zji Zj 2 B ajliajl „aj 2 ,a j2 ~ ( z j 1 B aji , aji + Zj 2 S« j2 ,» j2 ) + U , 

11a: ji ^ ki < j 2 = k 2 


c (k 1 , j 2 ', jl, J2) — Z kl Zj 2 i?a fcl ,a 31 ;aj 2 ,ai 2 Z kl £0^,0^ , 

lib: ji = ki < j 2 ^ k 2 

c \j 1) k 2 ', jl, J 2 ) = Zj x Zk 2 ,ajj; afc 2 , — z k 2 B ak ^^ aj ^ , 

He: ji < j 2 = k x <k 2 

^(j 2, k 2j ji, j 2 ) Zj 2 T ^k 2 ■^cxk 2 ,o,j 1 , 

lid : k 1 <k 2 = ji < j 2 

C\ki, ji, ji, j 2) Zj x z kl B a ^ a ^. ak ^ a ^ T z kl B akl>ot j 2 , 
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Ill: all ji, j r 2, k\, k 2 are different 


c(k i, k 2 , ji, j 2 ) = z k i z k2 B, 


’ ^^2 ’ ^32 ^ 


where 


B a , 0 -,a', 0 ' = jj {B at/3 B— B a ^> B a i :/3 ) , 


(B22) 


and each aq denotes the index of the chain containing the lepton line lj. 

The remained task is to find out a convenient expression of B a g. a i gi as a polynomial of 
w a . Expressing 1/ _1 U rs in Eq. (1B221) by Eq. (IB 1311 . B a ,p-,a',t 3 ' reduces to a polynomial 


B, 


a, 0 \ a !, 0 ' 


^ ^ ^ ^ (Carla's £as£a'r) ( ^0r'^0's' £/3s'£/3'r') 


l<r<s<n l<r'<s'<n 


x(-l) 


r+r'+sH-s' 


det (f/(rs|rV)) 


(B23) 


Here det ( U(rs\r's ')) is the determinant of the (n— 2) x (n—2) matrix {^/(rsIrV)^} obtained 
from s „ <n by eliminating the rth and sth columns and the r'th and s'th raws. 

If we define the function f r<s (x ) on the domain {1, • • • , n — 2} by 


fr<s(%) < 


X 

1 < x < r — 1, 

X + 1 

r < x < s — 2, 

x + 2 

s—l<x<n— 


(B24) 


the matrix elements U{rs\r's') xy are expressed in terms of U r » s » as 

U{rs\r s ) xy = Uf r<s ^j r , <s 0 y -). (B25) 

The application of the similar manipulation as that gives Eq. TOD to det ( U{rs\r's ')) yields 
det (t/(rs|rV)) = ^ 


a i<---<a „~2 


X ( ^ ] £ ( a ) Cai,f r<s (a(l)) ‘ ‘ ‘ £a n - 2 ,fr<s(<r(n-2)) 

l CGSn-2 


X 


£ ( <J ) £c*i,/ r , <s ,(<7(l)) ' ' ' ^a n -2,f r / <s i(cr(n-2)) I • 


l a - e &„-2 


(B26) 
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By inserting this expression into Eq. (IB23I) . is expressed in a simple form 


^ ^ ^a\ ' ' '1^a n —2 0(®j ® , Gi, j G n _ 2 ) @(/3, (3 , Gi, j G n — 2 ) j 

ai<-<a n _2 

(B27) 


with 


0(g, f3] Gi, • • • , G n _ 2 ) = E (-l) r+S+1 (Ca.rffta - £a,s£/3,r) 

l<r<s<n 

X ^ ] £ ( cr ) Coi,/r< s (o-(l)) ' ' ' Ca„_2,/r< s (o r (n-2)) • 


ff£S„_2 


(B28) 


The next task is to examine if @(g, (3\ Gi, • • • , G n _ 2 ) coincides with A{pt, /?, a±, • • • , G n _ 2 ). 
As was done to derive the first equality of Eq. 4Eni), we replace the sum over the permuta¬ 
tions of degree (n — 2) on the indices of loops with that over the permutations on the chain 
indices to get 

0(g, P; Gi, • ■ • , G n _ 2 ) = E (-l) r+S+1 {Ca,rCp,s ~ £a,s£/3,r) 

l<r<s<n 

X £ ( <J ) ^«a(l)>/r<s(l) ' ' ' £a<T(„- 2 )>/rO>(n- 2 ) 

cr£&n- 2 

= E E (-i) r+a+ M*) 

l<r<s<n cre ©„_2 

X ^c„(i), 1 ' ' ' £a CT ( r _i), r -1 £a, r £a CT ( r ), r+1 
X £a CT („_ 2 ),s-l S/3,s sa <T ( 8 _ 1 ),a+l ' ' ' £a CT („_ 2 ),n 

+ E E (-ir + - +2 £ (<r) 

l<r<s<n freS n _2 

X ^(i), 1 ' ' ' sa a ( r _i),r-l S/3,r ?a„ (r ),r+1 

X ^Oi a -( s _2) , a 1 SfV, 5 ^C^cr(s_ 1) i ^“(“l —2) ) ^ 

^ 1 £■{&) Ca,o-(l) £/3, cr(2) £«i,ct( 3) ' ' ' £a n _ 2 ,cr(n) 
cr£& n 

= A(a, (3, Gi, • • • , G n _ 2 ). (B29) 
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The third equality is given by using the fact that the sequences of permutations 

[chr(l)) j Cy<r(r—l)j Q^o-(r)) j Q!< 7 (s—2); / 3 , 1), j 0 ^a(n— 2)] 

I * /?, Cy,x(l )i 5 ^cr(r—1) ? ^cr(r) i i ®-cr(s— 2 )i ^<x(s—1)? ; ^<x(n—2)] 

| —* [cq /?, oq, • • • , a n -2] , 

[cy<r(l)) i QV(r—lp 3 i C^o-(r); j C^<r(s—2p Qq C^cr(s—1); j QV(n—2)] 

I * [tT, /?, CTo-(l), , C^(r(r—1); C^cr(r); j ®cr(s—2); QV(s—1); j QV(n—2)] 

t—> [a, / 3 , aq, • • • , a n _2] , (B 30 ) 

where the two steps in the first case produce signatures (—i) r + s + 1 and e(a) respectively, and 
the two steps in the second case produce (—l) r+s+2 and e(cr) respectively. They combine 
to form all possible permutations of degree n. In this way, we finally reach at a compact 
formula for B Q! p ]a > t p'] 

B a ,/3-,a',0' = ^ ] z ai ' ' ' z a n -2 

oi< ■<a„_2 

xA(a, a\ aq, • ■ • , a n - 2 ) A(( 3 , $, aq, ■ ■ ■ , a„_ 2 ) • (B 31 ) 


This formula provides a way to calculate B aj p. a i t pi algebraically. 

The cases (I), (IIa), (IIb), (He) (IId) contain -B q ,/3;«',/3'i where at least two of the chain 
indices coincide with each other. Actually such polynomials can be calculated much rapidly 
instead of using Eq. (HUD- From Eq. m for B aa and the definition (Hum B af) , B Qta]a ,,p 
can be written as 


B 0: a',/ 3 ' 


1 

u 


dU 

dw n 


^ ^ £a',rU(U )rs£/3',^ 


\r, s=l 


dlJ 

dw n 


^ ^ ^a',rU(U )r« u vU (U 
r , s, u , u=l 
n 

^ ^ £a',r(U )rs^/3' 




\r, s =1 

n 

U Ur(U - 1 ) 

r, s, w, v=l 

<9(17(7/-%) 


<97/ 




<9uy 


r, 5=1 


<9w n 


This yields an equality 


5, 


a, a ; a', (3' 



(B32) 

dB a t p, 

dw a 

(B33) 
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This equality also follows from the expression (IBSD for B a t pi and Eq. (1B3 ID . While Eq. (IB31I) 
requires the calculation of two A(a, a', cq, • • • , a n _ 2 ), the formula (IB33I) allows us to obtain 
B a a . a / p, simply by looking for the monomials containing w Q in B a i pi. The use of Eq. (IB33D 
is very efficient for the calculation of B a a . a i pi. Finally, these results can be translated into 
those for C tl through Eqs. W, dB2H>, and (1B221) . 
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